function [diff] = Equilibrium(Param0,varargin)

Cal = varargin{1}; CalValues = structvars(Cal);
for ctr0 = 1:size(CalValues,1); eval(CalValues(ctr0,:)); end

Param1 = exp(Param0)./(1+exp(Param0)).*High+(1-exp(Param0)./(1+exp(Param0))).*Low;

Par.alfa_A  = Param1(1)*10^2;
Par.alfa_B  = Param1(2)*10^2;
Par.gama_A  = 2;
Par.gama_B  = 2;
Par.alfa_eA = Param1(3)*10^3;
Par.alfa_eB = Param1(4)*10^3;
Par.gama_eA = 2;
Par.gama_eB = 2;
Par.phy     = Param1(5); 
Par.dlt     = 0; 
Par.lmd     = Param1(6)+1; 
Par.trf     = 5.5/100;
Par.trf_A   = Par.trf;
Par.trf_B   = Par.trf;
Par.kapa    = Param1(7); 

Par.bet     = 0.5; 
Par.eta     = 1-Par.bet; 
Par.psy     = 2;
Par.ro      = 0.01; 

Par.Q_m1    = 1;
Par.Q_0     = 1; 
Par.Q_1     = 1; 

Par.pi_dom = Par.bet*((1-Par.bet)/Par.eta)^((1-Par.bet)/Par.bet);
Par.pi_frn = Par.bet*((1-Par.bet)/(Par.eta*(1+Par.kapa)*(1+Par.trf_B)))^((1-Par.bet)/Par.bet);
Par.pi_frn_toA = Par.bet*((1-Par.bet)/(Par.eta*(1+Par.kapa)*(1+Par.trf_A)))^((1-Par.bet)/Par.bet);

Par.tau_A = tau_A;
Par.tau_B = tau_B;

Par.L_A   = L_A;
Par.L_B   = L_B;

Par.mbar = (length(Cal.y_init)-1)/2;
Par.mtot = 2*Par.mbar+1;
Par.madj = ceil((1-Par.bet)/Par.bet*log(1+Par.kapa)/log(Par.lmd))-1;
Par.madj = min(Par.madj,Par.mbar);
Par.mfactor = (1+Par.kapa)^((1-Par.bet)/Par.bet);

Par.ind_compet = 1;

Par.num_r   = 1/4;
Par.dt      = 1/32;
Par.burn_ss = 0;
Par.time_ss = (Par.burn_ss+250)/Par.dt;
Par.time0   = Par.burn_ss/Par.dt+0;
Par.time1   = (Par.burn_ss+(mid_year-init_year))/Par.dt;

% ==== Leapfrogging & Profit Functions ====

steps = (Par.mbar:-1:-Par.mbar)';

LF_aux = zeros(Par.mtot);
LF_pmf = zeros(Par.mtot);
for ctr_loc = 2:2*Par.mbar

	LF_aux = LF_aux+diag(ctr_loc*(diag(LF_aux,ctr_loc)+1),ctr_loc);

end
LF_aux = LF_aux+diag(1*(diag(LF_aux,1)+1),1);

lf_pmf = flip((1:2*Par.mbar)'.^-Par.phy/sum((1:2*Par.mbar).^-Par.phy));
for ctr_pmf = 2*Par.mbar:-1:1

    LF_pmf(1:ctr_pmf,ctr_pmf+1) = [lf_pmf(1:ctr_pmf-1); sum(lf_pmf(ctr_pmf:end))];
    
end

prft_fnc_A = [ones(Par.mbar-Par.madj,1)*(Par.pi_dom*L_A+Par.pi_frn*L_B);
              ones(2*Par.madj+1,1)*Par.pi_dom*L_A; zeros(Par.mbar-Par.madj,1)];
        
prft_fnc_B = [ones(Par.mbar-Par.madj,1)*(Par.pi_dom*L_B+Par.pi_frn*L_A);
              ones(2*Par.madj+1,1)*Par.pi_dom*L_B; zeros(Par.mbar-Par.madj,1)];
          
wage_dom   = [ones(1,Par.mbar+1+Par.madj) zeros(1,Par.mbar-Par.madj)];          

% ==== Auxiliary Matrices, Q transition ====

Trans.Qloc_idom = -1*eye(2*Par.mbar+1);
Trans.Qloc_idom(1) = Par.lmd-1;

Trans.Qloc_frn_R  = rot90(eye(2*Par.mbar+1));
Trans.Qloc_ifrn_L = rot90(-1*eye(2*Par.mbar+1));
Trans.Qloc_ifrn_L(end,1) = Par.lmd-1;

Qloc_edom = Trans.Qloc_idom;
Trans.Qloc_edom = Qloc_edom+LF_pmf.*Par.lmd.^LF_aux;
Trans.Qloc_idom = Trans.Qloc_edom;

Qloc_efrn_L = Trans.Qloc_ifrn_L;
Trans.Qloc_efrn_L = Qloc_efrn_L+flipud(LF_pmf);
Trans.Qloc_ifrn_L = Trans.Qloc_efrn_L;

Trans_innov = Trans.Qloc_idom+eye(Par.mtot);
Trans_fail  = rot90(Trans.Qloc_ifrn_L)+eye(Par.mtot);
Trans_exo   = -diag(steps~=0)+diag(steps(1:end-1)>0,1)+...
               Par.lmd*diag(steps(2:end)<0,-1);
           
Trans.Qtr_exo = Par.dlt*Trans_exo';

Step = {steps,Trans_innov,Trans_fail,Trans_exo,[prft_fnc_A prft_fnc_B],wage_dom};

% ==== Auxiliary Matrices, MU transition ====          
          
Trans.Muloc_idom = -1*eye(2*Par.mbar+1);
Trans.Muloc_idom(1) = 0;

Trans.Muloc_frn_R  = rot90(eye(2*Par.mbar+1));
Trans.Muloc_ifrn_L = rot90(-1*eye(2*Par.mbar+1));
Trans.Muloc_ifrn_L(end,1) = 0;

Muloc_edom = Trans.Muloc_idom;
Trans.Muloc_edom = Muloc_edom+LF_pmf;
Trans.Muloc_idom = Trans.Muloc_edom;
        
Muloc_efrn_L = Trans.Muloc_ifrn_L;
Trans.Muloc_efrn_L = Muloc_efrn_L+flipud(LF_pmf);
Trans.Muloc_ifrn_L = Trans.Muloc_efrn_L;

Trans.Mutr_exo = Par.dlt*(-diag(steps~=0)+diag(steps(1:end-1)>0,1)+...
                  diag(steps(2:end)<0,-1))';

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

% Steady State
% ============================================

% save Set_postestimation_use Par Step Trans y_init y_last % to be used in counterfactual exercises; save running the code at calibrated parameters

[rrsqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, mat_y, g_sqn, ind_neg, ...
            prft_fnc_sqn, wage_dom_sqn] = ...
                                    InterestPath_m(Par,Step,Trans,y_init);    

if isempty(ind_neg)
    
    disp(' ')
    disp(['          R in SS Problem - Returning !!! ' ...
        datestr(now, 'HH:MM:SS') ' !!!']); 
    disp(' ')
    
    diff = 10^5;
    
    return
    
elseif ind_neg == 1

    disp(' ')
    disp(['          R&D negative somewhere - Returning !!! ' ...
        datestr(now, 'HH:MM:SS') ' !!!']); 
    disp(' ')
    
    diff = 10^5;
    
    return
    
elseif ind_neg == 2

    disp(' ')
    disp(['          R transition negative or NaN somewhere - Returning !!! ' ...
        datestr(now, 'HH:MM:SS') ' !!!']); 
    disp(' ')
    
    diff = 10^5;
    
    return
    
end

rr_sqn = rrsqn(:,burn_ss+1:end);   

% keyboard

% 0) Fundamentals
% --------------------------------------------------------

prft_fnc_A = prft_fnc_sqn(1:size(prft_fnc_sqn,1)/2,:);
prft_fnc_B = prft_fnc_sqn(size(prft_fnc_sqn,1)/2+1:end,:);

wage_dom_A = wage_dom_sqn(1:size(wage_dom_sqn,1)/2,:);
wage_dom_B = wage_dom_sqn(size(wage_dom_sqn,1)/2+1:end,:);

% 1) Qualities
% --------------------------------------------------------

QA = Q_sqn(1:size(Q_sqn,1)/2,:);
QB = Q_sqn(size(Q_sqn,1)/2+1:end,:);

QwA = sum(wage_dom_A.*QA+flip(1-wage_dom_B).*QA*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet));
QwB = sum(wage_dom_B.*QB+flip(1-wage_dom_A).*QB*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet));

qq_sqn = [sum(QA)./QwA; sum(QB)./QwB];
                
wq_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(1,:).^(-Par.bet);
            (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(2,:).^(-Par.bet)];  

xk_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wq_sqn(1,:).^(1/Par.bet)*Par.L_B;
         (Par.pi_frn_toA/Par.bet)^(1/(1-Par.bet))./wq_sqn(2,:).^(1/Par.bet)*Par.L_A];

% 2) R&D and Values
% --------------------------------------------------------

X_iA = X_sqn(1:size(X_sqn,1)/2,:);
X_iB = X_sqn(size(X_sqn,1)/2+1:end,:);

V_iA = V_sqn(1:size(V_sqn,1)/2,:);
V_iB = V_sqn(size(V_sqn,1)/2+1:end,:);

X_eA = Xe_sqn(1:size(Xe_sqn,1)/2,:);
X_eB = Xe_sqn(size(Xe_sqn,1)/2+1:end,:);

V_eA = alfa_eA*X_eA.^gama_eA*(1-1/gama_eA);
V_eB = alfa_eB*X_eB.^gama_eB*(1-1/gama_eB);


% ===================== Moments ==========================
% ========================================================

% Growth 
% -------------------------------------------------------- 

IA = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet)-...
    sum(alfa_A/gama_A*X_iA.^gama_A.*QA)-...
    sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA)+...
    pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA)+...
    pi_frn_toA/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB)+...
    wq_sqn(1,:).*sum(QA)+...
    trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);  

IB = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet)-...
    sum(alfa_B/gama_B*X_iB.^gama_B.*QB)-...
    sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB)+...
    pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB)+...
    pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA)+...
    wq_sqn(2,:).*sum(QB)+...
    trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA);  

g_I = [(IA(2:end)-IA(1:end-1))./IA(1:end-1); ...
        (IB(2:end)-IB(1:end-1))./IB(1:end-1)];

g_IA = mean(g_I(1,1:time1));
g_IB = mean(g_I(2,1:time1));


% R&D Intensities (RD: level, RND = ratio)
% --------------------------------------------------------

RD_eA = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA);
RD_eB = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB);

RD_iA = sum(alfa_A/gama_A*X_iA.^gama_A.*QA);
RD_iB = sum(alfa_B/gama_B*X_iB.^gama_B.*QB);

% Exports - Imports (always and everywhere an IG phenomenon)
% --------------------------------------------------------

X_A = pi_frn*L_B*sum(flip(1-wage_dom_B).*QA)/bet.*wq_sqn(1,:).^(1-1/bet);
M_A = pi_frn_toA*L_A*sum(flip(1-wage_dom_A).*QB)/bet.*wq_sqn(2,:).^(1-1/bet);

X_B = pi_frn_toA*L_A*sum(flip(1-wage_dom_A).*QB)/bet.*wq_sqn(2,:).^(1-1/bet);
M_B = pi_frn*L_B*sum(flip(1-wage_dom_B).*QA)/bet.*wq_sqn(1,:).^(1-1/bet);
 
% Profits
% --------------------------------------------------------

PI_A = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet);
PI_B = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet);

% 6) Wages - Price Indeces (IG)
% --------------------------------------------------------

W_dA = pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA);
W_dB = pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB);

W_fA = pi_frn_toA/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB);
W_fB = pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA);

W_A = W_dA+W_fA;
W_B = W_dB+W_fB;

Wl_A = wq_sqn(1,:).*sum(QA); 
Wl_B = wq_sqn(2,:).*sum(QB); 
 
% 7) Production - GDP (Income Approach)
% --------------------------------------------------------

GDP_A = PI_A+W_A+Wl_A; % IA+RD_tA; % (= Expenditure Approach with NX = 0)
GDP_B = PI_B+W_B+Wl_B; % IB+RD_tB; % (= Expenditure Approach with NX = 0)


% 8) Ratios
% --------------------------------------------------------

RND_iA = RD_iA./GDP_A;
RND_iB = RD_iB./GDP_B;

RnD_A = mean(RND_iA(1:time1));
RnD_B = mean(RND_iB(1:time1));

XA_rat = X_A./GDP_A;


% ================ Objective Function ====================
% ========================================================

Targets = [mu_last; Growth_A; Growth_B; RnD_int_A; RnD_int_B; 
           Entry_A; EXtoY_A]; 

Moments = [mat_y(1:32,time1+1);
            [(1+g_IA)^(1/dt)-1 (1+g_IB)^(1/dt)-1]'; [RnD_A RnD_B]'; 
           mean(sum(X_eA(:,1:time1).*mat_y(:,1:time1))); 
           mean(XA_rat(1:time1+1))];

diff01 = (0.05-100*(g_IB-g_IA)*(1/dt))*1000*(100*(g_IB-g_IA)*(1/dt)<0.05);

diff_vec = Moments-Targets;

disp(' ')
disp(['          Outer Loop Finishing at ' datestr(now, 'HH:MM:SS') ' !!!']);
disp(' ')

weights = [ones(32,1)/32; 1; 1; ones(4,1)];  
diff = sum([(100*(diff_vec)).^2.*weights;
                diff01^2]);

% keyboard